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ABSTRACT 


The capabilities of an existing three-dimensional incompressible 
Navier-Stokes flow solver, INS 3D, are extended and improved to 
solve turbulent flows through the incorporation of zero- and two- 
equation turbulence models. The two-equation model equations are 
solved in their high Reynolds number form and utilize wall 
functions in the treatment of solid wall boundary conditions. The 
implicit approximate factorization scheme is modified to improve 
the stability of the two-equation solver. Applications to the 
three-dimensional viscous flow inside the 80 by 120 feet open 
return wind tunnel of the National Full Scale Aerodynamics 
Complex (NFAC) are discussed and described. 
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SECTION 1 : INTRODUCTION 


Three-dimensional, incompressible turbulent flows are frequently 
encountered in many industrially important flows, such as low- 
speed wind tvinnels or the space shuttle main engine. These flows 
represent a class of complex problems of practical importance. 

The purpose of this work is to offer an engineering tool which 
can take into account important viscous effects such as separated 
flow regions, secondary flows, and which can be flexible, 
accurate and computationally efficient. 

To fulfill this goal, an existing computer code, INS3D, 
developped by Kwak et al. (Ref. 1) for solving the 
incompressible, Navier-Stokes equations in a three-dimensional, 
curvilinear coordinate system, is applied. Before describing 
the applications, the features of the code are briefly 
outlined. 

The INS3D code makes use of the artificial compressibility 
concept, introduced by Chorin (Ref. 2). It consists in adding a 
time derivative of the pressure term to the mass conservation 
equation, and results in solving a system of hyperbolic 
equations. An implicit time-differencing procedure (Ref. 3) is 
adopted to march in time until a steady state is reached. 

The code has been applied to solve several laminar flows, (Refs. 
4, 5 and 6), and turbulence effects have been taken into account 
by means of a zero-equation model (Ref. 7) . In this approach, 
turbulence is modeled through the use of time averaging 
procedures which introduce unknown Reynolds stresses into the 
equations of motion. These stresses are then modeled by relating 
them to known mean flow quantities. 

In the present work, two eddy viscosity models are investigated 
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1) the Prandtl mixing length zero-equation model (Ref. 8) , and 2) 
the k-epsilon two-equation model (Ref. 9) . Eddy viscosity models 
have been extensivley studied for complex, two-dimensional flows 
(Refs. 10-12). Their extensions to complex three-dimensional 
flows are encouraging but deficiencies have been pointed (Ref. 13) 
in particular, for complex three-dimensional compressible flow 
fields with large separated zones. The assessment of eddy 
viscosity models in the INS3D code for complex three-dimensional 
incompressible flows remains to be done. The goal of this study 
is to make the INS 3D code operational for complex three- 
dimensional turbulent flows using eddy viscosity models. 

In the following sections, the equations of motion, the 
turbulence models, the numerical methods and wall function 
botindary conditions used in the simulations are discussed. These 
are followed by a representative application to the three- 
dimensional flow inside the 80x120 wind tunnel at NASA-Ames 
Research Center. 
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SECTION 2: EQUATIONS OF MOTION 


The basic equations used to analyze turbulent flows are the 
Reynolds-averaged incompressible Navier-Stokes equations. These 
equations expressed in Cartesian form may be written ; 



= 0 


(la) 


+ 

9u. u . 
11 

3p 

9t 

8x . 

~ Tx . 


3 

1 


9t . . 

(lb) 


According to Refs 2 and 3, equation (la) is modified as follows; 


^ +8^-0 da') 

3t 3Xj_ 

where t is time ; are the Cartesian coordinates, u^^ are 
corresponding velocity components, p is the pressure and is 

the viscous stress tensor and 3 is the pseudo compressibility 
parameter. 


The viscous stress tensor includes both molecular and Reynolds- 
averaged turbulent contributions. By means of the eddy-viscosity 
hypothesis, the viscous stress tensor is written in the following 
form; 


ID 


= -(V + 


Vt> 


(u 


i/D 




( 2 ) 


where v and , are the molecular and turbulent (eddy) 

viscosities, and svibscript notation has been introduced for 
partial derivatives. 
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As mentionned in the introduction, the eddy viscosity for the 
two-equation models is represented in terms of additional field 
variables here taken to be s^ and s^ so that = v^(s^,S 2 )* 
The additional field equations governing s^^ and s^ are written 
as 


^®i 3 

3^ + ^ <®i"j ^ ‘^ij) 


«i 


qii = - (V + 

^ri 


(3) 

i = 1,2 
j = 1,2,3 


with = (X/y,z), 

where the modeling constants, Pr^^ , Pr 2 are turbulent Prandtl 
numbers, and , H 2 are source functions. The source functions 
depend on s, , s,, and u. . , and will be described in detail in 
the next section. 


The equations are written in a dimensionless form with: 



k= U^ic 


X = Lx ~ UUj V = ULv 

^ (4) 



For applications to curvilinear coordinates {^,r\,0 , equations (1) 
are transformed by means of the following relations : 


C = Ux,y,z) 

n = n(x,y,z) (5) 

C = C(x,y,z) 

We will also use the condensed form : = ?,ri,C i =1, 2, 

or 3. The tildes C") have been dropped for convenience. The reference 
length is L, the reference velocity is U. 
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Further details on this transformation may be found in Ref. 1. 

The governing equations in conse3rvation-law form are expressed in 
generalized curvilinear coordinates as : 


/N 

9t 


+ 


' where 



Q = 3/j = -i 


- = 0 
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_ — 
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u,- 
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'^l 

'' 1 
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"2°i + ^i2P 
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U 3 U. + L. 3 P 

• « 


— — 


1,2,3 


( 6 ) 


(7) 


=vi = 1 


1,U2,U33 


J is the Jacobian of the transformation and 


°i = ^il“l + ^i2“2 + ^i3“3 




= (5i)^ ; 


l^ 2 - (Ci)y ; Li3 - (5i>, 


are the contravariant velocities and the metrics of the 
transformation respectively. 

Equation (7) is simplified to : v . ) 9 T (8) 

J ^ ^ 

The approximation (8) is valid for nearly orthogonal grids. The use 
of such approximation results from a compromise in the choice of 
saving computer time and storage. In the application presented 
here, the choice of keeping available the maximum nimiber of grid 
points was made. Therefore, the metrics are not kept in memory 
and are computed when they are needed and the approximation (8) 
is made to minimize computing time. 
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SECTION 3: TURBULENCE MODELS 


In this section we will describe the turbulent models used in the 
present study. As mentioned in the Introduction, these are the 
zero- and two-equation eddy viscosity models. Advantages and 
disadvantages of eddy-viscosity models have been studied (Refs. 11 
12, 13 and 14) , some of these aspects are recalled here. 

3.1 Zero-Equation Model 

A relatively large number of zero-equation models have been used 
with the Navier-Stokes equations. These are the simplest of all 
turbulent models because they require no additional field 
equation and contain only a few modelling constants. The 
nvunerical compatibility of these models is high, which is why 
they have been widely used. However, zero-equation models are 
deficient in two important aspects: First, the assumption on the 
algebraic length-scale that they require is very difficult to 
specify for complex separated flows, which reduces the generality 
of the model. Secondly, these models do not account for flow 
history or stress relaxation effects, which are important in 
complex flows. However, these models have given reasonable 
agreement for simple flows, and they are useful for comparison 
purposes in evaluating and improving performances of more 
sophisticated models. Here the simplest and the oldest zero- 
equation model is employed, i.e. the Prandtl mixing length model. 
This is a two layer model with a modified distance function. An 
earlier investigation using a modified distance function has been 
performed by Hving and MacCormack (Ref. 15). The eddy viscosity 
and related parameters are written below: 

= 1^ |Vxu| (9) 

The eddy viscosity is related to the local absolute value of 
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vorticity, which is invariant with respect to a rotation of the 
coordinate system, and involves a single unknown parameter: the 
mixing length 1. The specification of this parameter is given 
below. 

The Prandtl mixing length, is given by : 

1 = min ( K 1^ , 1^ ) with <=0.4 (10) 


The inner length or distance function is 1^^ 


d d^ 
n ? 




where d^ is the distance from the wall ^ = 0, and d^ is the 
distance from the two (side and top) walls, corresponding to n=0. 


The definition of the inner length is critical for the evaluation 
of the turbulence length-scale. This composite formula is designed 
to account for the size of turbulence eddies or the turbulence 
mixing length near the corner under the influence of both walls. 
Note that : 

If 

If i.^d^ 

The outer length is 1 q = 0.09 6(x) , where 6(x) is the 
boxindary layer thickness. A simplified formula for 6(x) is used 
in the present study and is obtained from the empirical flat 
plate formula (Ref. 16) : 


6 (X) ^ 0.37 ; 

X 


R 


ex 


Ux 

V 


X is the streamwise distance from the inflow boundary. 

This simplified formula for 6 (x) has been used in the present 
study because of the great difficulty in computing the boundary 
layer thickness in the conventional manner by locating the edge 
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of the boundary layer. Since the main purpose of this work was to 
implement and make operational the two-equation model in the 
INS3D code, this simplification was used as a short cut to reach 
this goal in a short time. 
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3.2 Two-Equation Model 


In these models, two field equations are specified (see Eq. 3) . One 
equation determines the velocity scale, the other the length-scale. 

In this work, the high Reynolds number form of the k-e. model is 
used in conjuntion with wall function boundary conditions. The 
variable k represents the turbulent kinetic energy, the variable t 
is the dissipation rate of kinetic energy. This model is 
summarized below, further details of this model can be found in 
Ref. 10. 

(k-e) model; s^ = k S 2 = e 


= P - . 

H, = Cl p . Cj 


Making use of the gradient diffusion hypothesis, the term 
describing the production of turbulence, P, is written as : 


P = [vJu. ■ _ 

T 1,3 3,1 

The constants in the model are : 




Cy = 0.09 


= 1.45 


P 


rk 


1.0 


P 


re 


1.3. 


C 2 = 1.92 


The choice of these constants are discussed in Reference 17. 
Briefly, Cy is determined from the observation that the 
production and the dissipation in a constant stress layer flow 
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adjacent to a wall are balanced, convection transport is 

negligible. The constant C2 is obtained from grid turbulence 

results where diffusion and production are zero. Knowing these 

two constants and the von Karman constant, can be calculated 

from the e equation which is greatly simplified in near wall 

region. The diffusion constants or Schmidt numbers, P and P 

rk re 

were determined from systematic comparisons of computations with 
experiments for a series of flows. 
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SECTION 4; NUMERICAL METHODS 


The nvimerical method used to solve the incompressible Navier- 
Stokes equations is detailed in Ref. 1 , it is an implicit, 
approximately factored finite-difference method due to Beam and 
Warming (Ref. 18) . The same numerical procedure to solve the two- 
equation turbulence model is employed and is described in Ref. 10. 

An important feature of the numerical method not considered in 
Ref. 10 is the treatment of the turbulence field equations, in 
particular, the implicit treatement of the turbulence source 
functions. This treatment follows the work by T.J. Coakley (Ref. 

14) and is given here; such a treatment is helpful to the 
successful numerical implementation of two-equation models. In 
the present approach, the continuity and momentum equations are 
solved first and the primitive variables p, u, v, w are updated 
before going on to solving for the two turbulence variables. 

Therefore, p, u, v, w are known at the time level n+1 when the 
resolution of the two-equation model starts. Let us consider the 
two-equation model system of partial differential equations 
expressed in curvilinear coordinates : 

D^ + E-+F +G =S (13) 

t 5 n c 

D = (k/J, e/J) ; E, F, G are the turbulent fluxes associated with the 
directions 4, C. This is written as: 

AD n+1 

^35 9n 3c ■ ^ 

In the INS3D code, continuity and momentum equations are solved 
first with the pressure and velocity components updated before 
solving for the turbulent variables k and e . Therefore, 
the turbulent source terms are expressed at time n+1 instead of 
at time n like in the basic algorithm. Then, Taylor expansions 
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give ; E 

F 

G 

S 


A = 


3E 


,n+l 

1 

= 

E^ + 

AE = 

,n+l 

— 

F^ + 

B AD 

,n+l 

r 

= 

G^ + 

C^AD 

,n+l 

> 

= 

+ 

K^AD 

1 

1 ' 

B 

3F 

3D 

/ c 


3D 


n n 
e" + A AD 


3G 

3D 


3S 

3D 


As a result, a Jacobian Matrix K, derived from a Taylor expansion 
of the source term appears in the Left Hand Side (LHS) of the 
finite-difference equations. We write : 

[1 + At( 3-A + 3 B + 3^C - K)] ad = -At C 3 ^E^ + 3 + 3 - S^] 

5 n ? V § n ^ 

Rf?S ' 

The term K is a rather complicated Jacobian matrix which will be 
approximated by the identity matrix multiplied by a constant. 

The unfactored algorithm is : 

Cl - At K + At (L^ 

where L^, L^, are three one-dimensional operators. The 

algorithm is expressed in this convenient form to use the spatial 
splitting technique, characteristic of the alternating direction- 
implicit (ADI) scheme. Then, we write 


(I - At K) Cl + At (1 - At K)"^ (L^ + L + L^)] AD = RHS 

? n c 


The operators L^, L^, become : 


= (1 - At K)"^ 


etc. . 


Then, the approximate factorization ADI method is used, and the 
implicit source term involving the constant becomes product of 
the time step; 


is factored as ; 


(1 + At L^) (1 + At L^) (1 + At L^) AD =-At(l - At K)"^RHS 
The choice of K follows -the choice of Ref. 14. We set: 


At— 


At 


1 + At K 

A composite formula is used for K : 
K = a (2C2 § ) + (1 - a) 2 ^t ) 


1_ ^ 
AS^ 


AS^ AS^ AS^ 

5 n ? 


,.2 2 . 
+ 


y\ + 


AS' 


and 


n n -^n n 


s. = 


2 


2 

Yr + 


'C ' 


X, = 2L(J+1.)..- x(J-i) ^ 
" 2 


For the linear problem, this method is unconditionally stable. 

This algorithm has not been fully analyzed in the present work. 

For example, various combinations of the constant a could be 
tried to determine its optimtun value, (i.e. a =0, 1., 0.5 
etc..). In this work, the value of a that performs well is a = 0. 
The advantage of this algorithm over the case with K = 0, is the 
enhancement of the stability. 
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SECTION 5; WALL FUNCTION BOUNDARY CONDITIONS 


The use of wall functions was originally advocated by Launder and 
Spalding (Ref. 19) and applications have recently been extended 
to complex compressible flows (Ref. 20) . 

The simplest boundary conditions to solve the Reynolds-averaged 
Navier-Stokes equations are no slip conditions. However, when 
turbulence models are used, this simplicity is challenged by the 
need of special treatment of low- Reynolds niimber terms, to solve 
the viscous sublayer. This treatment, called the integration to 
the wall, varies accordingly to the turbulence model in use (see 
Ref. 20) . The alternative approach, used here, is the use of 
"wall functions", or the use of the logarithmic law of the wall 
formula for the velocity profile near the surface. The reasons 
for this approach are driven by a concern of economy and improved 
stability of the turbulence model resolution. The use of wall 
functions eliminates a considerable number of mesh points , and 
the need to capture the rapid changes in mean flow and turbulence 
quantities in the near wall region. Another reason of this choice 
is that the use of wall functions has been demonstrated to 
perform well for complex flows (Ref. 20) , in particular, the 
predictions of separated flows with the law of the wall are in 
some cases improved over the same predictions using the 
integration to the wall approach. 


The high Reynolds number form of the two-equation model is 
employed in conjuntion with the law of the wall procedure, 
similar to the one developed and tested in the TURF code (Ref. 
11) for compressible flows. This formulation is adapted to the 
INS3D code as follows : the friction velocity u^ is detemined 
implicitly from the empirical formula (14), given Y 2 ' ^ 


u. 


^2 " 


In 


EUry2 


E = 9.128 


(14) 


K V 

tho subscript 2 denotes the first point from the wall. 


21 



. 2 

The wall shear stress t is obtained from the friction velocity: x„=u^ 

w w T 

We assume : ^ * 

Letting = Yj =0 (no slip) enables to compute the effective eddy 
viscosity at the mid point 1% , which gives 


The molecular viscosity v is sometimes neglected in this formula. The 
eddy viscosity mid way between the wall and the first point from 
the wall is needed to evaluate the viscous fluxes. Two ways were 
employed : 


1. This eddy viscosity is obtained using an average viscosity : 


^2 ^w 
^TUs “ U2 

- V =-i- 
2 

*“t1 '"t2 

2. Another approach is to assume 

: = 2 

(^2Im - v) 


"2 


^T2 

The first procedure demonstrated to be the best approach, the friction 
velocity reaches its final value faster. 


The eddy viscosity v ^2 computed from the basic formulas for 
each model: 


2 -»■ 

a) zero-equation model: v ^2 ^ 

b) two-equation model: v _2 = ~ 


In the case of the zero-equation model, one sided differencing in the 
direction normal to the walls are used to compute the vorticity at the 
first mesh point off the wall. 


In the case of the two-equation model, boxindaty conditions on the two 
turbulence variables are needed, these are discussed page 32 of this 
report . 
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SECTION 6: RESULTS 


The resulting nximerical model was applied to the simulation of 
the flow inside the 80 by 120 feet wind tunnel located at NASA 
Ames Research Center. A general schematic of the wind tunnel is 
shown in Figure 1. The design of this wind tunnel is the result 
of a major effort in aerodynamic design. The theoretical part of 
this activity was accomplished with the use of singularity-type 
panel computer code (Refs. 21 and 22) . A three-dimensional Euler 
analysis of this flow has been performed by Kaul et al. (Ref. 

23) , showing results in good agreement with the experimental 
data. In the present study, viscous effects are taken into 
account through the use of turbulence models. 

Three computations are reported and compared: 

1. A laminar computation with a Reynolds niimber of 1000. 

2. A turbulent computation using the Prandtl mixing length model 
and the law of the wall boundary condition, at a Reynolds number 
of 6.6x10^ . 

3. A turbulent computation with the two equation Jones Launder 

model and the law of the wall boundary condition, at a Reynolds 

7 

number of 6.6x10 . 


6 . 1 Grid Generation 

The computational grid shown in Figures 2 through 8 is the result 
of intensive efforts performed through trial and error. 

The three-dimensional grid is a combination of algebraic 
procedures and available elliptic and parabolic grid generation 
computer programs. Figure 2 represents the ground plane, the 
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inlet part was generated by the parabolic grid generation program 
of T.A. Edwards (Ref. 24) . The nozzle part was generated by an 
elliptic grid generation program (Kristin Hessenius ' private 
commianication) . Figure 3 is an enlargement of the inlet region in 
the ground plane, with the K-lines emerging perpendicularly from 
the side wall (K = 1) . This two-dimensional grid is rotated 
around the K=KMAX axis to discretize the full computational 
domain. Figure 4 shows the grid in the vertical plane (L = LMAX- 
1) . Figure 4 a is an enlargement of the inlet lip area, generated 
by the parabolic grid generation program (Ref. 24) . Figure 5 is 
the cross section grid downstream in the rectangular test 
section. Figures 6 a and b represent two perspective views of the 
three-dimensional grid, (a) viewed from the wind tunnel entrance 
and (b) viewed from the exit. The vertical plane, L = LMAX-1, is 
a plane of symmetry. The advantages of this grid are : 

1. It can be mapped into a cube. 

2. Points can be clustered only close to rigid boundaries. 

3. A zonal approach is not necessary. This grid results in 

a one block grid and avoids the complications of patching grids. 

4. Orthogonality of the grid lines is accomplished at the rigid 
boundaries . 

A disadvantage is that the line K = KMAX is a singularity line 
where all the L- lines are merging together (at the bottom wall in 
the symmetry plane) . This singularity caused some difficulties in 
applying numerical boundary conditions, and this will be reported 
later. 

This grid was used for the laminar computation. It employs 93 
planes in the streamwise direction (JMAX = 93) , 45 planes in the 
spanwise direction (KMAX = 45) , and 25 rotated planes (LMAX = 25) . 
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For the turbulent computations, the grid points had to be 
redistributed to satisfy the criterion of validity of the law-of- 
the wall. When applying the law of the wall boundary conditions, 
special attention must be given to the first spacing from the 
wall, since the law of the wall is valid for : 

20 < y"*" < 1000 approximately 

where y+ “t ^ “x^ex'^ ! R = 22- 

” “ 5 — = “55 '' 

d is the distance from the wall of the first mesh point, x is the 
streamwise distance from the grid line J = 1. 

For a turbulent flow above a flat plate: u Pf t 

"Iff " (Ref. 16) 

Therefore, the spacing of the first mesh point is defined by the 
algebraic formula : 

d = 5.8 X R“* V 
ex 

After trial and error , the first spacing of the final 

4. 

"turbulent" grid is obtained with an initial y of 600, resulting 
in a first spacing from the wall test section of 0.001 approximately. 

The K and L grid lines were redistributed as shown in Figures 7 
and 8, the final grid size is ; 

JMAX = 83, KMAX = 33, LMAX = 33. 


6.2 Laminar solution 

The goal of the present study is to perform turbulent simulations 
of the 80 by 120 feet wind tunnel, however, it is necessary to 
begin with a laminar simulation that will help defining initial 
flow conditions, boundary conditions, time steps and that will be 
used as a comparison with the turbulent computations. 

The Reynolds number is 1000 and is based on half the width of 
the wind txinnel (equal to 60 feet) . 
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Initial conditions. 


The dimensionless pressvire is set to 1 everywhere in the flow 
field. The initial values of the velocity are scaled to conserve 
the mass flux rate at the outflow. 

Boundary conditions. 

At the walls: zero-gradient extrapolation for the pressure and no- 
slip condition for the velocity. 

At (singularity point) K=KMAX: no slip boiandary conditions for 
velocity components and pressure were employed at first. Then , 
from comparisons of computations with experiments, it appears that 
two point extrapolations of the velocity components were the most 
appropriate boxindary conditions. 

At the inflow boundary: the inflow boundary is placed sufficiently 
far upstream so that the velocity gradients are small here. The 
stagnation pressure and flow direction (taken parallel to the grid 
line directions) are prescribed. 

At the outflow: a simple one point extrapolation is used for the 
velocity components and pressure. This is the most stable 
boxindary condition. However, because this boundary condition does 
not automatically conserve the mass, at the start of the 
computation mass weighting procedures on the pressure and 
velocity, as described in Ref. 4 are employed. 

The results of this computation are shown in Figures 9 through 
16. The velocity vector plots in the streamwise direction are 
shown in Figure 9 for the plane L=16 (which runs diagonally accross 
the test section) . A separation bubble immediatly downstream of 
the cowl, was observed experimentally but is not predicted by the 


26 


computation. This is probably attributed to a too coarse grid 
spacing. The choice of a coarse spacing in the inlet area was 
made so that sufficient grid points could be used in the area of 
the test section. Figure 10 shows the crossflow velocity vectors 
at X = 210 ft, at the beginning of the test section. Figure 11 
shows the crossflow velocity vectors at x= 240 ft. 

Figure 12 shows the variation of the axial velocity (normalized by 
the test-section velocity) with the spanwise coordinate, y, at an 
X station situated 27 ft upstream of the inlet, close to the 
z = 0 plane. Indicated on this Figure are the results of the 
inviscid analysis (Ref. 23) , referred to as "inviscid" on this 
Figure and the following ones. Experimental measurements are 
indicated by circles. We observe that the inviscid analysis seems 
to agree better than our viscous analysis, however, this location 
X = -27 feet, has the coarsest grid spacing, coarser than the 
inviscid analysis. Also, the viscous results are obtained by 
three-dimensional interpolation which introduces an absolute 
error of about 5%. Figure 13 shows a variation of the normalized 
axial velocity with the vertical coordinate z at x = -27 ft, near 
the symmetry plane between the two side walls of the tvinnel. 

Figure 13 shows the INS 3D results in disagreement with the 
experimental data and the inviscid solution, at the time of this 
computation, this was attributed to the use of a constant eddy 
viscosity model (laminar) that is known to overpredict the 
thickness of boundary layers. However, the same discrepancy was 
observed with the turbulent solutions, action to correct it 
was taken at that time and this will be described in the next 
subsection. We will see that this is caused by the use of no-slip 
condition at K=KMAX, and the use of two point extrapolation 
formula for the velocity components resulted in improved 
agreement of the computation with the experimental measurements. 

Figure 14 shows resultant velocity along the center line (y=z=0) 
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versus downstream distance from the cowl to the test section. The 
results are compared with the panel method and with the inviscid 
method (Ref. 23) . The inviscid theory results are for the 
velocity at the ground. The laminar results are for the velocity 
mid way between the groiand plane and the top plane, so that they 
can be compared with inviscid theories. In this case, the laminar 
velocity variation agrees with the panel method results better than 
with the U. Haul's results Figure 14 also shows similar comparisons 
for the tiinnel floor pressure distribution versus the downstream 
distance from the cowl to the test section. We observe some 
discrepancy in the inlet cowl area. Downstream the inviscid and 
laminar results are in better agreement. 

Figure 15 presents the pressure variation with the downstream 
distance on the side wall centerline near the z = 0 plane. The 
comparison between the computations and the experiments is good 
except near the separation bubble where the pressures given by 
the laminar simulations are underpredicted with respect to the 
experiments and inviscid solutions. 

Discussion: Figure 11 has shown the existence of weak 
longitudinal vortices in the test section, one close to the floor 
center-line, the other close to the top surface. The origin of 
this vorticity in the flow was not clear and a concern was that 
this vorticity could be created at the inflow by a poor inflow 
boundary condition. Figure 16 shows how they originate 
by presenting velocity vectors in the boundary layer at L=3, L=5 
and at the edge of the boundary layer and it is concluded that 
these vortices result from viscous effects. Because of the side 
wall contraction, a negative pressure gradient is created 
crosswise, this deflects the fluid toward the center, the slow 
moving fluid of the boundary layer is deflected faster than the 
fluid in the main "inviscid" flow. Downstream of the contraction, 
the streamline curvature changes and the situation is reversed. 
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Figure 16 a shows the downstream velocity vectors in the boundary 
layer directed toward the side wall. Figure 16 c shows the 
downstream velocity vectors outside the boundary layer directed 
toward the center line, this effect causes the flow to rotate. 
These vortices have been studied experimentally (Refs. 25-27) , 
they will occur in laminar boundary layers as well as turbulent. 
Other than that, the tendency of the flow to be directed toward 
the center of the plane is observed, as this was also seen in 
flow-visualization studies. 

6 . 3 Zero-equation model solution 

The results obtained with the zero-equation model solution are 
shown in Figures 17 through 27. 

The zero-equation model computations were initialised from the 
same solution than the one employed for the laminar computation. 
The boundary conditions were at first unchanged from the laminar 
solution. 

Figures 17 through 21 present quantitative results compared with 
inviscid theories and experiments. Figures 22 through 27 present 
qualitative results as vector plots. 

Figure 17 presents the normalised axial velocity variation with 
the spanwise distance at x = -27 feet. We observe good agreement 
with the experimental data, close to the wall and in the center 
of the wind tunnel (y=0) . Mid way the INS3D results overpredict 
the axial velocity variation, however, it should be recalled that 
the grid spacing is particularly coarse in this region and this 
could explain this apparent discrepancy. 

Figure 18 presents the normalized axial velocity variation with 
height z at x = -27 feet. The dotted line indicates the results 
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obtained when a no~slip boundary condition was used at K=KMAX, 
resulting in a disagreement with the inviscid solution and 
experiments. At first, it was believed that the inflow boundary 
condition was the cause of this problem, but various inflow 
velocity profiles were tried (with a boundary layer and without a 
boundary layer) , resulting in the same results. Then, the 
boundary condition at K=KMAX was modified instead of using no- 
slip condition, zero-gradient in the z-direction (actually K 
direction) was employed, i.e. all quantities at K=KMAX for all L 
planes are set equal to the corresponding value at K=KMAX-1 and 
L=LMAX-1. This provided improved agreement, the INS3D results are 
the solid line in Figure 18. From now on, all the results 
presented were obtained with this boundary condition. 

Figure 19 presents the resultant velocity variation with the 
downstream distance, the INS3D velocity results are evaluated in 
the center of the wind tunnel on the axis of symmetry so that 
they can be compared with the inviscid theories. Figure 20 shows 
the tunnel floor centerline -pressure variation with the 
downstream distance. Indicated in Figures 19 and 20 are the 
results of the inviscid analysis (Ref. 23) and of the panel 
method. We observe some discrepancy of the pressure variation in 
the inlet lip area, while the three computed results agree 
downstream of the inlet tunnel. 

Figure 21 shows the pressure variation on the side wall 
centerline with the downstream distance compared with the 
inviscid analysis and experiments. Here too, we observe some 
discrepancy of the pressure variation in the inlet lip area and 
this is attributed to the deficiency of the computation in 
not predicting the separated flow in this area. We think that 
this is caused by a too coarse mesh and the zero-equation model 
itself. 


Figures 22 and 23 present results of velocity vector plots in the 
plane L=22 and L=15. Figure 24 is an enlargement of the velocity 
vector plot in the plane L=15 near the inlet cowl, we observe 
that the velocity profiles in the most curved part of the inlet 
cowl seem to be near separation. 

Figures 25 a, b, c show crossflow velocity vector plots at three 
locations, x=32 feet, x= 164 ft, x= 228 ft. Figures 25 a and b 
show the flow converging toward a point on the symmetry plane. 
Figure 25 c indicates the existence of a weak vortex in the 
middle of the test section. As a comparison, we chose to show 
Figure 26: crossflow velocity vectors at the same location than 
the one of Figure 25 c. The results of Figure 26 were obtained 
with the no-slip botindary condition for the velocity at K=KMAX 
and we observe that the qualitative features of the flow are 
quite different from the one observed on figure 25 c. The 
computation using no-slip boundary condition for the velocity at 
K = KMAX predicts two vortices, one close to the floor, the other 
close to the top, while the computation using zero-gradient 
extrapolation for the velocity indicates the presence of a single 
vortex, in the center of the computational domain. 

Figure 27 shows the velocity profile along the line K=22, in the 
center of the test section. 

6 . 4 Two-equation model solution 

Figures 28 to 38 shows the results obtained with the two-equation 
model . 

The same computational grid as the one used for the zero- 
equation model solution was utilized. The two-equation model 
solution was restarted from the zero-equation model solution, 
with free stream values of k and e used for initial conditions. 
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Free stream values. An important aspect of the k and solution 
is the determination of the free stream values, and this is 
discussed here. Free stream values of the two turbulent variables 
appear to control the stability of the solution. If the free 
stream values are too large, the turbulence kinetic energy and 
the turbulence dissipation cannot diffuse, the eddy viscosity 
remains null and the solution is laminar. If these values are 
too small, and in particular, the value of the turbulence 
dissipation, the eddy viscosity becomes vinrealistically large, 
causing the computation to become unstable. These values were 
determined by trials and errors and finally, the ratio of eddy 
viscosity to laminar viscosity was set equal to 0 . 1 , with 
k = 7 l 0 xl 0 ”^, e is deduced from these two values and from the 
Kolmogorov formula. 


Boundary conditions. 

The boundary conditions on the pressure and velocity components 

were unchanged from the one used for the zero-equation model 

solution. Associated with the use of wall functions, boundary 

conditions on the two turbulence variables are needed. Two 

alternative wall boundary conditions on k and e were tested: 

2 
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y is the distance from the wall 


2 . k 
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and 


» k2 = 


K 


(extrapolation from interior) . 


, nd <16) 

£9 = eoYo/yo = , 3 indicates the 2 mesh point 

from the wall. 


These boundary conditions are critical in getting a good start of 
the computations. Attempts to start computations with the 
boundary conditions (16) failed. At the start of a computation, 
the friction velocity is small, the k and e values are small 
also, and with the boundary conditions (15) , k and e profiles 
develop smoothly, however, it was observed that the friction 


32 


velocity obtained from the law of the wall does not grow as fast 
as the k and e values do. As a result, a small discontinuity in 
the profiles of k and e was observed and although the computation 
remains stable, the accuracy of the k and e quantities near the 
wall was in question. At that point, a switch to boundary 
conditions (16) was made to insure a smooth monotonic growth of 
the k and e profiles, and the friction velocity is allowed to 
adjust to its steady flow value. After about 1000 cycles, the 
solution is close to convergence, and boundary conditions (15) 
and (16) provide similar results. 

While performing these computations, a serious stability problem 
was encountered and is disccused here. It was found that, once 
the k and e profiles have developed, the solution could not 
remain stable. The problem is caused by e reaching its free 
stream value too quickly at the fifth mesh point from the wall, 
while the k profile retains a monotonically decreasing behaviour. 
This results in an eddy viscosity too large. Many attempts to 
stabilize the solution failed. These attempts consisted in trying 
various combinations of the boundary conditions on k, e and on 
the eddy viscosity, and in using a small time step. At that 
point, the stability of the k-e solution algorithm was in 
question, and the following test was performed. The two-equation 
model solution was restarted with the converged zero-equation 
model solution. The eddy viscosity solution given by the zero- 
equation model was "frozen" while all the other variables, 
pressure, velocity and turbulence variables were computed. It 
was observed that the k and e profiles develop smoothly, the 
computation remains stable, and a relatively large time step 
could be used. From this test, we deduced that the k and e 
algorithm was numerically compatible and that our problem was an 
inaccuracy on the calculation of e caused by a too coarse mesh. 

A test was introduced in the computation ; if e is less or equal 
to its freestream value and k is greater than its freestream 
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value, e was set equal to : 




£ 


K 




This formula is based on the experimental measurements of 
Klebanoff for flat plates (Ref. 28) . This was sufficient to 
provide e with a smooth transition from its computed value to a 
fixed freestream value. The effect was rather spectacular and 
stabilized the computation. However, this is only a temporary fix 
and given a finer grid spacing, this condition would probably 
become unnecessary. 

Figures 28 through 32 are a set of quantitative results for the 
two-equation model, similar to the ones shown for the zero- 
equation model. We observe that the two-equation model solution 
provides agreement with inviscid theories and experiments similar 
to the results obtained with the zero-equation model. In 
particular, in Figure 32, we observe that in the inlet cowl area 
the pressure variation is vinderpredicted by the two-equation 
model and exhibit a relatively strong variation. To improve 
agreement in this area, we believe a finer grid spacing would be 
helpful. 

Figures 33 through 38 show velocity vector plots illustrating the 
qualitative features of the flow field. Velocity vectors in the 
plane L = 22 are shown in Figure 33. Figure 34 presents an 
enlargement of the velocity vector plot in the plane L=22 near 
the inlet cowl, the predicted boundary layer thickness is thin, 
no separation is observed and this agrees with experimental 
observation. Velocity vectors in the plane L = 15 are shown in 
Figure 35. Figure 36 is an enlargement of the velocity vector 
plot in the plane L=15 near the inlet cowl where separation was 
observed experimentally. Here, we observe a stronger gradient 
variation in the boundary layer in spite of the coarse mesh 
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(velocity vectors for all grid points are plotted ) , the velocity 
profile in the cowl center appears to be close to separation. 

Given a finer spacing, separation might be predicted, as observed 
in experiment. 

Figures 37 a and b show velocity vector plots in the cross- 
section of the wind tunnel at x = 32 feet and x = 228 feet. At 
X =228 feet, we observe a weak vortex in the middle of the test 
section, caused by the lateral divergence of the streamlines 
which intensifies the lateral component of vorticity. 

Figure 38 compares eddy viscosities and velocity profiles 

obtained with the zero- and two-equation model at x= 210 feet 

along the diagonal line L=22. We note here that the boundary- 

layer thickness predicted by the two-equation model was somewhat larger 

than the one given by the zero-equation model. With only 10 

grid points in the boundary layer, the coarse mesh inhibits an 

accurate prediction of the boundary layer thickness. For 

instance, at x = 210 feet along the line L = 22, the boundary 

layer thickness is such 


7.2 feet < 6 

< 

9.6 feet 

(zero-equation model) 

9.6 feet < 6 

< 

12 feet 

(two-equation model) 


Figure 39 shows the turbulence kinetic energy and turbulence 
dissipation curves at x = 210 feet, along the line L=22. In this 
Figure, a rapid decrease in the dissipation toward the outer edge 
of the boundary layer is observed, this sharp decrease occurs 
within three mesh points from the wall, causing the eddy 
viscosity to become large. Experimental observation (Ref. 28) 
confirms this behaviour. 
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SECTION 7: COMPUTER EXECUTION TIME 


The computations reported here require approximately: 

a) laminar case: 0.0001 seconds per grid point per iteration, 

b) zero-equation model case: 0.00011 seconds per grid point per 
iteration, 

c) two-equation model case: 0.00014 seconds per grid point per 
iteration. 

Typically, 700-1200 iterations are required for laminar 
solutions, and about 1500 iterations are required for turbulent 
solutions. However, the two-equation model solutions reported 
here, were run up to 3000 iterations to gain confidence in the 
stability of the computation but 1000 to 1500 iterations appear 
sufficient to reach a converged solution. In the two-equation 
model solution, the residuals of the two turbulence variables 
decrease quickly (within 500 cycles) by three order of magnitude 
and then decrease slowly by an approximate factor of 0.0002 per 
cycle. 

For the original algorithm, the largest time step that maintains 
the stability of the code was 0.003. The time step used to start 
computations was 0.0001, and was increased progressively to 0.003 
up to 700 cycles. The zero-equation model and the two-equation 
model solution employs the same time step. The laminar solution 
was obtained with a coarser grid and the largest time step was 
0.007. The algorithm modification described in Section 4 was not 
fully investigated because of lack of time, however, it was 
observed to perform best after about 500 cycles, (the time step 
could be increased faster). At the start of a computation no 
significant improvement on the size of the time step was noted. 

We feel that this algorithm modification is promising and more 
in-depth investigation needs to be done. 
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SECTION 7; CONCLUDING REMARKS 


This section sxBiimarizes our conclusions, the first subsection 
outlines the accomplishements of this effort, the second 
subsection attempts to state objectively the deficiencies of this 
numerical simulation and gives recommendations for future 
improvements of viscous three-dimensional simulations. 

7 . 1 Accomplishements . 

1. Zero- and two-equation turbulence models have been 
incorporated into the INS3D code, including special 
boundary conditions appropriate to the use of wall 
functions. This procedure enables the use of coarser 
mesh spacing near the wall, compared with integration to 
the wall procedures (where no slip boundary conditions are 
used) and leads to greater stability and efficiency. 

2. Laminar and turbulent solutions (with a zero-equation 
model and with a two-equation model) were computed and 
the results were compared with available experimental 
measurements and with inviscid solutions. 

3. These computations simulate a flow in which viscous 
effects, as measured by extent of separation, were small 
so that the results obtained using the INS 3D code should 
be (roughly) comparable to the results obtained using 
inviscid codes. This is what was observed. Relatively 
minor differences from inviscid results were noted for 
the pressures and axial velocities. 

4. The modifications to improve the stability of the 
algorithm was introduced , but not extensively tested. 
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For the present, the results are inconclusive regarding 
improvements in convergence rate. 

5. A new simple-single block three-dimensional grid was 
developed. However, some improvements could be made, in 
particular, by increasing the number of mesh points. 


7 . 2 Recommendations . 

1. The numerical algorithm modification for improving the 
stability of the computations was not optimized, and more 
work (following the outlines given in Section 4) is 
required to make this modification fully efficient. 

2. The solutions in the inlet region do not predict the 
separation observed in experiments, this is probably due 
to a lack of resolution. A finer grid could cure this 
problem, or a new type of grid might be required. The 
best possible grid is an orthogonal single block grid and 
a three-dimensional elliptic grid generator code would be 
very helpful. 

3. The grid singularity at K=KMAX caused some difficulties 
regarding the use of boundary conditions. Although 
changing from no slip to zero-gradient type of boundary 
conditions imroved the solutions , more improvement is 
needed, because of the retardation of velocity profiles 
as the singularity is approached. This could also be 
possibly caused by the mesh discontinuity along the 
diagonal or by the use of the thin layer Navier- Stokes 
equations rather than the full Navier-Stokes equations. 

We think that it is essential to use the full Navier- 
Stokes equations in the predictions of three-dimensional 
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viscous flows, associated with a more regular grid. 


4. The two-equation model solutions are not entirely 

satisfactory, a special stabilizing fix was needed to 
achieve stable solutions. Here too, more work is needed 
on wall funtion boundary conditions to remove this problem. 

We established the importance of using different models (constant 
eddy viscosity (laminar) , zero-equation and two-equation models) 
combined with experimental comparisons. Such a process is 
essential to provide reliable predictions, in Computational Fluid 
Dynamics. 
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(365 ft) 



REF = 60 feet. 



- 200.0 - 175.0 - 150.0 - 125.0 - 100.0 - 75.0 - 50.0 - 25.0 
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Figure 4. Planar grid in the y-plane 
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Figure 5. Cross-section plane grid, x = 210 feet. 
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Figure 6 a. Grid plot showing the surfaces J=1 , L=1 , L=LMAX. 
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Figure 6 b. Grid plot showing the surfaces 



=KMAX 
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Figure 9. Velocity vectors in the plane L =16. 
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Figure 10. Velocity vectors in the cross-section plane, x= 210 feet. 
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12. Normalized axial velocity variation with the spanwise 
distance y at x = -27 feet. 
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Normalized axial velocity variation with height z, 
at X = -27 feet. 
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14. Centerline velocity (resultant) and tunnel floor 

centerline pressure variation with the downstream distance. 


56 





15. Pressure variation on the side wall centerline with 
the downstream distance x. 
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16 a. Velocity vectors in the olane L 
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16 b. Velocity vectors in the plane L = 5. 
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17. Normalized axial velocity variation with the spanwise 
distance at x = -27 feet, zero-equation model. 
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O PANEL METHOD AT Y = 5 ft 
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Velocity resultant variation with the downstream 
distance, zero-equation model. 




Tunnel floor centerline pressure variation with the 
downstrecun distance, zero-equation model. 
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Pressure variation on the side wall centerline with 
the downstream distance x, zero-equation model. 
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22. Velocity vectors in the plane L = 22, zero-equation 
model. 
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near the inlet cowl, zero-equation model. 
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25 a. Crossflow velocity vectors at x = 32 ft, zero-equation 
model. 
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25 b. Crossflow velocity vectors at x = 164 ft, zero-equation 
model. 
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25 c. Crossflow velocity vectors at x = 228 ft, zero-equation 
model. 



Crossflow velocity vectors at x = 228 ft using no-slip 
boundary conditions at K = KMAX. 
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28. Normalized axial velocity variation with the spanwise 
distance y, two-equation model. 
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29. Normalized axial velocity variation with height z, 
two-equation model. 
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30. Velocity resultant variation with the downstream 
distance, two-equation model. 
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31. Tunnel floor centerline pressure variation with the 
downstream distance, two-equation model. 
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32. Pressure variation on the side wall centerline with 
the downstream distance x, two-equation model. 


78 


33. 


P^OjC!/ i3 
OF POOR QUALITY 


A 


A 


\ 


-4 


Re = 6.60 X 10^. 

RRin = R*? Y RR Y 



’‘/Lref 


Velocity vectors in the plane L 


22, two-equation model. 
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37 a. Crossflow velocity vectors at x = 32 ft, two-equation 
model. 
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Comparisons of eddy viscosities and velocity profiles for 
the zero- and two-equation models. 
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